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Collimated Jet or Expanding Outflow: Possible Origins of GRBs 

and X-Ray Flashes 
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o : 

^ i ABSTRACT 

3 : 

We investigate the dynamics of an injected outflow propagating in a progeni- 
■ tor in the context of the collapsar model for gamma-ray bursts (GRBs) through 

two dimensional axisymmetric relativistic hydrodynamic simulations. Initially, 
>- ■ we locally inject an outflow near the center of a progenitor. We calculate 25 

models, in total, by fixing its total input energy to be 10 51 ergs s -1 and radius 
of the injected outflow to be 7 x 10 7 cm while varying its bulk Lorentz factor, 
q ■ r = 1.05 ~ 5, and its specific internal energy, eo/c 2 = 0.1 ~ 30 (with c being 

speed of light). The injected outflow propagates in the progenitor and drives a 
large-scale outflow or jet. We find a smooth but dramatic transition from a colli- 
CLi 1 mated jet to an expanding outflow among calculated models. The opening angle 

O 1 of the outflow (6> sim ) is sensitive to T ; we find 9 sim < 2° for T > 3. The maxi- 

^ ■ mum Lorentz factor is, on the other hand, sensitive to both of T and eo; roughly 

Tmax ~ r (l + eo/c 2 )- I n particular, a very high Lorentz factor of r max > 100 is 
achieved in one model. A variety of opening angles can arise by changing eo, even 
^ ■ when the maximum Lorentz factor is fixed. The jet structure totally depends on 

r . When r is high, a strong bow shock appears and generates a back flow. 
High pressure progenitor gas heated by the bow shock collimates the outflow to 
form a narrow, relativistic jet. A number of internal oblique shocks within the 
jet are generated by the presence of the back flow and/or shear instability. When 
r is low, on the contrary, the outflow expands soon after the injection, since the 
bow shock is weak and thus the pressure of the progenitor gas is not high enough 
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to confine the flow. Our finding will explain a smooth transition between the 
GRBs, X-ray rich GRBs (XRRs) and X-ray Flashes (XRFs) by the same model 
but with different e values. 

Subject headings: hydrodynamics - jet - GRBs - supernovae - shock - relativity 

1. INTRODUCTION 

The GRBs are to our best knowledge the most energetic phenomena in the Universe. So 
far intense efforts have been made both on the observational and theoretical grounds toward 
understanding of their natures, but their origins still remains to be investigated. One of the 
most important finding recently is that GRBs are involved with relativistic collimated flows. 
To account for the observations, an extremely high bulk Lorentz factor, typically more than 
100, is required (Rees & Meszaros 1992). 

GRBs are known to be composed of two classes: long-duration GRBs (with duration 
being longer than a few seconds) and short- duration GRBs (with duration being less than a 
second). At least, some of the long-duration GRBs are known to be associated with super- 
novae (SNe). Good examples of GRB-SN connection are GRB980425/SN1998bw (Galama 
et al. 1998) and GRB030329/SN2003dh (Stanek et al. 2003; Uemura et al. 2003; Hjorth et al. 
2003; Price et al. 2003). These provide strong evidences that the central engines of (at least 
part of) the long-duration GRBs are SNe. Such association was theoretically predicted by 
Woosley (1993) and Paczynski (1998). A signature of a supernova contribution is actually 
found in the afterglow spectra of these GRBs. These supernovae are categorized in the type 
Ic whose progenitor has lost its hydrogen and helium envelope when the core-collapse occurs. 
Another possible GRB associated with SN is GRB021211. A strong absorption feature is 
observed in the spectrum of the afterglow (Delia Valle et al. 2003). They concluded that 
the absorption is due to Call which is synthesized by the associated supernova explosion. 
Further, the discovery of the host galaxy of the long duration GRBs being star forming 
galaxies (Bloom, Kulkarni, Djorgovski 2002; Le Floc'h et al. 2003; Tanvir et al. 2004) also 
strengthens the idea of strong GRB-SN connection. 

It has been suggested that even supernovae, in which no associated GRB was found, 
might have a link with GRB. For example, the peculiar SN2002ap recorded a high velocity 
component of 0.23c and huge kinetic energy of the jet of 5 x 10 50 ergs, at least. These values 
are similar to those of GRBs, indicating a similar explosion mechanism of SN2002ap to that 
of GRBs. Totani (2003) concluded that SN2002ap is one example of the supernovae which 
failed to make a GRB. 
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Recently, a number of X-ray-rich GRBs (XRRs) and XRFs, very similar phenomena 
to GRBs but with significantly lower peak energy, have been successively discovered thanks 
to the good performance of HETE-2 (see Heise et al. (2001) and Sakamoto et al. (2005)). 
Interestingly, the event rates of XRRs and XRFs are similar to that of the long duration 
GRBs (Sakamoto et al. 2005). The origin of these events is poorly understood, but similar 
burst properties of GRB, XRRs, and XRFs except for peak energy leads to an idea that they 
all might have the same origins but with different viewing angles (Nakamura 2000) or with 
variable opening angles (Lamb et al. 2005). 

In this paper, we elucidate the theory of XRFs, XRRs, and GRBs in the context of 
the so-called collapsar model. The collapsar is a death of a massive star in the last stage 
of the stellar evolution. The collapsar model for GRBs was proposed by Woosley (1993; 
see also MacFadyen & Woosley 1999), for a central engine of GRBs. In this model, strong 
outflows, or jets, emerge from deeply inside the collapsar and propagate into the interstellar 
medium (ISM), producing gamma-ray bursts. MacFadyen and Woosley (1999) performed 
two dimensional hydrodynamic simulations based on this model. Assuming annihilation of 
neutrino and anti-neutrino, they deposited thermal energy around the center of the core 
which had been collapsed and become a black hole. Their initial mass density profile is 
very flattened due to rotation of the progenitor. The gas around the center, where the high 
thermal energy is deposited, expands and forms very collimated outflow; i.e., a "jet". The 
outflow successfully became a bipolar outflow. Unfortunately, however, their calculations 
were not relativistic one, so the relativistic effects which are important to understand GRBs 
were not included. 

Since the mass density of the progenitor is quite high, it is not a trivial issue whether 
or not the formed outflow can always keep collimated structure and break out from the 
progenitor surface as a jet. It has been pointed out through the AGN jet simulations that the 
multi-dimensional effects are so important for the dynamics of the "light jet" into some dense 
gas (see, e.g., Mizuta et al. (2004), and references therein). Here, light jet stands for the jet, 
the mass density of which is smaller than that of the ambient gas. At least two dimensional 
hydrodynamic calculations are indispensable to investigate the outflow propagation and its 
dynamics inside and outside the progenitor. 

Multi-dimensional, relativistic hydrodynamic simulations have been so far performed 
by several groups in the context of collapsar model (Aloy et al. 2000; Zhang, Woosley, & 
MacFadyen 2003; Zhang et al. 2004; Umeda et al. 2005). Aloy et al. (2000), for example, 
performed relativistic hydrodynamic simulations based on the model by MacFadyen and 
Woosley (1999). They have found that a bipolar flow is created and it breaks out from the 
progenitor. Interestingly, the maximum Lorentz factor of about 40 has been achieved, when 
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the jet breaks out of the progenitor. Another type of relativistic hydrodynamic simulations of 
the outflow have also been performed by Zhang, Woosley, & MacFadyen (2003), Zhang et al. 
(2004), and Umeda et al. (2005). They modeled a mainly very hot jet, i.e. an initially thermal 
energy dominated jet. Injected jets from the computational boundary, which is assumed to 
be very close to the collapsed center of the progenitor, always successfully propagate in the 
progenitor keeping good collimation and break out the progenitor. 

It is still open question, however, whether a collimated outflow emerging from the center 
of the progenitor can always break out or not. Although several provenance studies have 
demonstrated successful propagation throughout the progenitor and breakout of the input 
outflow, they might have assumed unrealistically large energy input in the initial condition. 
We should be aware that the formation mechanism of the outflow from the center of the 
collapsed progenitor has been poorly understood. In other words, we still do not know 
the physical conditions (density, thermal energy, kinetic energy, magnetic energy, opening 
angle etc) for generating outflows. Further, it is not completely clear yet what discriminates 
between SNe associated with GRBs and those without GRB association. The connection 
between XRFs and GRBs is another important issue. A key factor may be attributed to the 
different dynamics of the outflow propagation. 

Motivated by these questions we perform series of relativistic hydrodynamic simulations 
of the outflow propagation in the progenitor and ISM. In these simulations, we fix the total 
input energy power of 10 51 ergs s" 1 but vary the bulk Lorentz factor (r ) and the specific 
internal energy of the initial outflow (e which excludes the rest mass energy). We discuss 
what types of outflow can emerge from the central system for a wide range of parameters, 
r and eo- 

This paper is organized as follows. In Sec. 2 we introduce our model and then explain 
the numerical methods and the initial background and outflow conditions. The results of the 
numerical simulations are presented in Sec. 3, where we will demonstrate the emergence of 
two distinct types of outflows: a collimated jet and an expanding outflow. We then discuss 
the dynamics and structure of the outflow, focusing on the distinctions between the two 
types of flows in Sec. 4. The final section is devoted to conclusions. 

2. NUMERICAL METHOD AND CONDITIONS 

2.1. Our Model: Progenitor and Injected Outflow 

Our model is based on the collapsar model (see, e.g., Woosley 1993). According to 
this model the release of the gravitational energy is the main energy source of the energetic 
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phenomena. This model explains the formation of an outflow as follows: When an iron core 
of the massive star collapses, a system consisting of a black hole or a proto-neutron star and 
an accretion disk is formed at the center of the progenitor. Some fraction of the collapsing 
gas can produce an outflow, although the formation mechanism of this outflow is not fully 
understood yet, as in the case of AGN jets. Annihilation of neutrinos and anti-neutrinos 
emitted from the accretion disk or via MHD process is a possible scenario to generate such 
outflows. The progenitor is expected to spin rapidly. The gas along the rotational axis can 
freely fall into the center while the gas along the equatorial plane only gradually falls because 
of stronger centrifugal force. As a result, a tenuous regions are created along the rotational 
axis. After the core collapse, the outer envelopes begin to fall freely. The free fall timescale 
of the envelopes is longer than the dynamical time scale for the jet to propagate within the 
progenitor and to hit its surface. After the formation of the outflow, it should propagate in 
the progenitor and break out the surface of the progenitor into ISM. Finally the outflow is 
observed GRB and an afterglow. 

We adopt the radial mass profile from Hashimoto (1995), assuming that the progenitor is 
spherically symmetric when the iron core was collapsed (see however Petrovic et al. (2005); 
Yoon & Langer (2005); Woosley & Heger (2006) for recent calculations on the effects of 
the angular momentum and magnetic fields). At least near the core, the profile along the 
equatorial axis is very similar to that used by MacFadyen and Woosley (1999) and Aloy et 
al. (2000). The progenitor had a mass of about 40 solar masses in the main sequence and has 
16 solar masses in the pre-supernovae stage. The hydrogen envelope has already been lost 
by the stellar wind. Although our progenitor includes helium envelope which would produce 
type lb SN, main properties of the outflow dynamics, such as, collimated jet or expanding 
outflow described below, strongly depend on the mass profile around the injection point, i.e. 
silicon, carbon, and oxygen envelopes. Our main conclusion is not affected whether we adopt 
the progenitor which includes helium envelope or not. 

Figure 1 shows the radial mass density profile from the center to the surface of the 
progenitor. The mass density decreases from 10 10 g cm -3 to 1 g cm -3 from the center to the 
surface located at 3.7 xl0 10 cm from the center. The lower boundary of the computational 
domain (z\ ow ) is set to be 2 x 10 8 cm from the center of the progenitor. The center of the 
progenitor is set to be z = (see Figure 2 for a schematic view of the progenitor and the 
location of the computational box). Note that the total mass of the progenitor within this 
radius is about two solar masses; that is, here we postulate the situation that material of 
about two solar masses collapses towards the center and forms the system of a proto-neutron 
star or black hole and a surrounding accretion disk. This distance between the center of the 
progenitor and the computational lower boundary is the same as that adopted by Zhang, 
Woosley, & MacFadyen (2003). We will also calculate an additional case study, in which the 
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lower boundary is set to be at 2 x 10 7 cm from the center, following Aloy et al. (2000), in 
which thermal energy is deposited as a driver of the outflow. 

Since the pressure becomes very high due to the strong bow shock, the pressure of the 
progenitor can be negligible. The initial thermal energy in the progenitor is set to be very 
low. 

The formation mechanism of the outflows around the core is not fully understood yet. 
We assume injection of an initial outflow from the lower boundary and that the direction 
of the initial outflow is parallel to that of the cylindrical (z) axis. This method is basically 
the same as that adopted by Zhang, Woosley, & MacFadyen (2003) Zhang et al. (2004) 
and Umeda et al. (2005), although the injected outflows calculated by Zhang, Woosley, 
& MacFadyen (2003) and Umeda et al. (2005) had a finite opening angles and were not 
initially parallel flow. In the present study the radius {Ro) and power (E ) of the injected 
outflow is fixed to be Ro = 7 x 10 7 cm and E = 10 51 ergs s -1 , respectively. In this paper 
the subscript stands for the values of the injected outflow and E does not include the 
rest mass energy. The energy flux (E /irR 2 ) does not depend on radius. The net injected 
energy during the first ten seconds amounts to 10 52 ergs. This values is close to the explosive 
energy of SN1998bw and SN2003dh and is higher than the normal explosive energy of ~ 10 51 
ergs (Iwamoto et al. 1998; Woosley, Eastman, & Schmidt 1999; Hjorth et al. 2003). Such 
energetic supernovae are sometimes categorized as hypernovae. 

Two more parameters are necessary to specify the outflow condition. The bulk Lorentz 
factor (To) and specific internal energy (eo) are chosen in this paper. These parameters 
characterize the kinetic and thermal energy per particle in the injected outflow, respectively. 
As eo and/or To increases, so does the kinetic and/or thermal energy per particle. The rest 
mass density and pressure of the injected outflow can be derived from these two parameters, 
assuming an equation of state. The larger eo and/or T is, the smaller becomes the rest mass 
density. For example, the rest mass density (po) is given as 

p c 2 = E [((l + e /c 2 7)r 2 - ro)^]" 1 ^]" 1 , (1) 

where 7 is adiabatic index. The ideal gas equation of state, p = (7 — l)pe, is employed, 
where p is pressure, and the definition of specific enthalpy h/c 2 (= 1 + e/c 2 + p/pc 2 ) is used. 
The equation (1) is derived from T 01 — poro^o = E/{ttRq). 

Figure 3 shows rest mass density in the (T , e ) plane for fixed total energy, E = 
10 51 ergs s" 1 and the radius, R = 7 x 10 7 cm, of the injected outflow. The symbols in the 
plane present the calculated models. 

Assuming that all the initial total energy is efficiently converted to kinetic energy, we 
can estimate the maximum bulk Lorentz factor from the energy conservation law, E tot 
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M b c 2 r max (Piran 1999) where M b is baryon mass of the outflow. Approximately the total 
energy is the sum of the rest mass, kinetic, and thermal energy, E tot ~ M^c 2 [l + (r — 1) + 
r e /c 2 ]. Then the maximum Lorentz factor of a fire-ball can be estimated as (Piran 1999) 

r m ax~r (i + eo/c 2 ), (2) 

which means almost all thermal energy is converted to kinetic energy during expansion like 
SN explosion (Arnnet 1996). 

We will systematically vary T and e to survey any possible physical situations, since 
we are unaware which is the case. The condition which determines these parameters should 
depend on the formation mechanism of the outflow from the center of the progenitor, i.e., 
we should know in detail how the mass and the angular momentum are distributed in the 
progenitor, how a massive star collapses, how the neutrino transport occurs, and how the 
magnetic field grows and affects the dynamics. For example, a rapidly rotating massive star 
can form geometrically thick accretion disk around a new-born black hole or proto-neutron 
star after the core-collapse. The amount of the neutrino and anti-neutrino emission from such 
an accretion disk could be large. Then the annihilation rate of neutrino and anti-neutrino 
will be large enough so that thermal dominated outflow may be generated. If the accretion 
disk is thin due to a lack of rapid rotation of the progenitor, conversely, the annihilation rate 
will be small, which may results in the formation of a baryon rich outflow. Note, however, 
that these conclusions may not be firm, since we still do not know the formation mechanism 
of this outflow. Magnetic fields could be a key here. 

Since there is a possibility that a non-relativistic but collimated outflow is generated from 
the center of the progenitor, we explore both the relativistic and non-relativistic outflows. 
The bulk Lorentz factor, r , is varied from the relativistic to the non-relativistic regimes, 
r = 5, 4, 3, 2, 1.4, 1.25, 1.15, 1.1, and 1.05 (which correspond to v /c ~ 0.98, 0.97, 0.94, 0.87, 
0.7, 0.6, 0.5, 0.4, and 0.3, respectively). The specific internal energy, e /c 2 , is changed from 
30, 5, 1.0, and 0.5. As a result, the rest mass density of the outflow is in the range from 
2.4g cm" 3 [the fastest and hottest outflow : (r ,e /c 2 ) = (5.0,30)] to 4 x 10 4 g cm" 3 [the 
slowest and coldest outflow : (t> /c, eo/c 2 ) = (0.3,0.5)]. The pressure of the outflow also 
changes. Table 1 summarizes 25 calculated models. 

To compare with Zhang, Woosley, & MacFadyen (2003) we also list the ratio of the 
total energy to kinetic energy, f , is presented in the table, where the rest mass energy is 
excluded from the total energy; 

y o = p ro(r - 1) ^ 

Po(h /c 2 )T 2 - (Po/c 2 ) - p r ' 
(see Zhang, Woosley, & MacFadyen (2003)). Because the definition of f in Zhang et al. 
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(2004) was changed to be the ratio of kinetic energy to total energy, we also show the inverse 
of /o in the table. Our model A300, (r , eo)=(5,30) is similar one studied by Zhang et al. 
(2004); Umeda et al. (2005). The theoretically estimated Lorentz factor T max from Eq. (2) 
is also shown in Table 1. Some cases, for example (vq/c, eo/c 2 ) = (0.3, 5.0), produce subsonic 
flow and are, hence, difficult to calculate. We exclude such cases from Table 1 to assure 
good numerical accuracy. Since the mass density of the progenitor in the innermost region 
for the computation is about 10 6 g cm" 3 , all our models produce initially so-called light jet. 
If the flow can keep the collimated structure in the progenitor, such an outflow is expected 
to interact with the back flow (Mizuta et al. 2004) and has complex internal structures in 
the jet. From eq. (2) we understand that the most predominant case for the GRBs is the 
outflow with very high T and high e ; i.e., model A300 (r max ~ 150). The second one is 
model A50 (r max ~ 30). 



2.2. Hydrodynamic Equations 

We numerically solve two dimensional relativistic hydrodynamic equations, assuming 
the axisymmetric geometry. The basic equations are, 

djpT) + 1 dr{pYv r ) + d{ P Yv z ) 

dt r dr dz 
d(phT 2 v r ) 1 dr(phT 2 v 2 + p) d(phT 2 v r v z ) 

dt r dr dz 

d(phT 2 v z ) 1 dr(phT 2 v r v z ) d(phT 2 v 2 + p) 

dt r dr dz 

d(phT 2 -p) | 1 dr(phT 2 v r ) | d(phT 2 v z ) 

dt r dr dz 

V{ (i = r, z) is the i-th velocity component. The equations are written in the unit that the 
speed of light is unity. The updated version of numerical hydrodynamic code used in Mizuta 
et al. (2004) is employed in this study. The code adopts Godunov-type scheme which is 
advantageous for capturing strong shocks with a few grid points in good accuracy. In this 
version the physical values, such as, pressure, rest mass density and the space components 
of 4-velocity, are used for this interpolation using the Piecewise Parabolic Method (PPM) 
which allows us to get third order accuracy in space. The version of third order accuracy in 
space is used. The results of ID and 2D test calculations are presented in the Appendix. 

In this study, we calculate the propagation of the outflow crossing through the progenitor 
for 10 seconds after the outflow injection. Since the outflow crossing timescale is shorter than 
the free-fall timescale, we neglect the gravitational field by the core. We also ignore self- 
gravity, for simplicity. When an outflow propagates from the inner boundary towards the 
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(5) 
(6) 
(7) 
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surface of the progenitor, a strong bow shock appears and drives the progenitor gas to high 
pressure and temperature. Nucleosynthesis could occur in the such a hot and high density 
medium. Since produced entropy due to the nucleosynthesis is much smaller than that by 
the strong shock, we do not consider nucleosynthesis as an energy source in the energy 
equation (Eq. 7). We assume the ideal gas equation of state, that is p = (7 — l)pe, where 
7(= 4/3 constant in this study) is adiabatic index. 

We used non-uniform grid points. Logarithmically uniform 500 grid points are spaced 
for 2 x 10 8 cm < z < 6.6 x 10 10 cm. We also set uniform 120 zones for < r < 1.2 x 10 9 cm 
and logarithmically uniform 130 zones for 1.2 x 10 9 < r < 1.1 x 10 10 cm. Some of our results, 
especially those having slower injection velocity (vo < 0.7c), exhibit spreading outflow, which 
has a large transverse component of the flow with respect to the direction of the injected flow. 
We find in the test calculations of ID shock tube problem with relativistic transverse velocity 
(see, A. 2) that the limited resolution of our numerical codes produces some numerical errors 
(numerical errors in both positional and absolute values of density, pressure, 3-velocity, up 
to several tens of percent). However, in some slower models, such as, models G01, H01, 
and 101, the flow is close to non- relativistic flows and the errors caused by the relativistic 
transverse velocity are negligible. In some faster models, such as, models A300, A50, A10, 
A05, etc, the relativistic flow in mainly appear in the collimated jet and propagates along 
the cylindrical axis. Since the transverse velocity is not so large, the numerical errors caused 
by the transverse velocity are also negligible. We can thus conclude that the errors in mildly 
relativistic and spreading models may affect to the opening angle, but does not alter our 
main conclusion of a smooth transition from the collimated jet to the expanding outflow 
that is discussed in Sec. 3. The achieved maximum Lorentz factors in each simulations are 
not affected by this numerical error. 

The boundary condition at z = 2 x 10 8 cm is reflective except the inner 7 grid points 
which are used to inlet outflow. The reflective boundary condition is also employed at 
cylindrical axis. The outflow boundary condition is employed at r = 1.2 x 10 9 cm and 
z = 6.6 x 10 10 cm. We compellingly replace the numerical flux at the boundary using 
injection conditions, such as, po, €q, and T , to inlet the outflow. The numerical fluxes of the 
mass, momentum for z direction, and energy at the injection points (z = z\ ow , < r < Rq) 
are given as p ro^o; Po^orV^o 2 + Po, an d poh T 2 v , respectively. 
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3. COLLIMATED JETS AND EXPANDING OUTFLOW 

3.1. Collimated Jets: Cases with High r 

First, we show the results of the high Lorentz-factor cases with T = 5; i.e., models A01, 
A10, A50, and A300. Figure 4 shows the contours of rest mass density (upper) and Lorentz 
factor (lower) in each model at the time of breakout of the jet; t = 3.50 s (A300), t = 3.50 s 
(A50), t = 3.00 s (A30), and * = 2.75 s (A01), respectively. Note that the lower-left corner of 
each panel does not correspond to the center of the progenitor (see again Fig. 2). The region 
with a large Lorentz factor is localized and is found only along the cylindrical axis. This 
means that the outflow propagates throughout the progenitor, keeping a very collimated 
structure. Thus we call this collimated outflow a jet. The half opening angle (6> sim ) of the 
jet in each model is kept small; it is only 9 sim ~ 1°, when the breakout occurs (see, Table 1). 

The collimated jet is surrounded by the shocked progenitor gas and back flow. The 
width of the jet is about 10 9 cm when the jet breaks. Lateral shocks are found at distances 
of about 5 — 7x 10 9 cm from the z-axis. Since the mass density of the injected outflow is much 
smaller than that of progenitor near the center (see the previous section), the velocity of the 
head of the jet is less than that of the jet (~ c), as was analyzed by Norman, Winkler, & 
Smarr (1983) and Marti et al. (1997). The jet takes about 3 sec to cross over the progenitor. 

The collimated jet eventually breaks out of the surface of the progenitor. We observed 
small differences in the crossing times of the jet within the progenitor among these models. 
The jet in model A01 passes through the progenitor in a shorter time than that in model 
A10, and, the jet in model A10 passes in an even shorter time than that in model A50 
and A300. This can be understood, since a lighter jet, i.e., lower density jet, feels larger 
resistance to proceed in the progenitor than a denser jet. 

The dynamics and morphology of the calculated jet, such as its collimated structure 
and the appearance of a back flow, are very similar to those shown by Zhang et al. (2003, 
2004). The flow structure is also very similar to those simulated in the context of AGN jet 
propagation, although most previous simulations have been done under the assumption of a 
constant ambient density. One of the most prominent features of the jets in high-r models 
is the appearance of the back flow (see Fig. 5a). In this figure, only main flow, i.e., a jet, 
and back flow from the head of the jet are shown. 

Since the jet propagates in the progenitor whose gas density decreases outward, the 
hot spot does not clearly appear, to the contrary to the simulation of the jet propagation 
into a dense gas cloud with no density gradient. At the head of the jet three discontinuities 
exist: bow shock or forward shock (FS) which drives progenitor gas to a high pressure and 
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high temperature state, contact discontinuity (CD), and reverse shock (RS) or terminal 
Mach shock. The kinetic energy of the outflow is converted to thermal energy thorough the 
reverse shock. The gas which passed this reverse shock forms a back flow in the anti-parallel 
direction to that of the jet. The path of the back flow is not straight but is meandering. 
The interaction between the main jet flow and the back flow enhances the oblique shocks in 
the jet flow. As a result, the profile of the physical quantities, such as density, pressure, and 
Lorentz factor, is not monotonically decreasing nor increasing outward. Figure 6 shows one 
dimensional profile (along the z axis) of the rest mass density and Lorentz factor of models 
A300, A50, A10, and A01 in the early phase of the simulation. The radial pressure profiles 
exhibit frequent up and down, like a teeth of a saw, so does the density profiles except at 
around the head of the flow. Some discontinuities are seen along the z-axis, although these 
features are difficult to see in the contours in Fig. 4. In model A10, A50, and A300 the 
Lorentz factor also shows a teeth-of-saw structure. Those discontinuities correspond to the 
internal oblique shocks. 

Before the breakout the bulk Lorentz factor increases during the propagation in the 
progenitor only up to about 43(A300), 20 (A50), 9.6 (A10), and 5.6 (A01), respectively. 
After the breakout, the bulk Lorentz factor further grows up 41 in model A50 and 104 in 
model A300 (see Fig. 7 for a series of one dimensional profiles of the rest mass density and the 
Lorentz factor along the cylindrical axis of model A300). The high Lorentz factor (r > 100) 
originates from the components which are injected in the later phase (t > 7 s), although the 
head of the outflow has already passed through the boundary (z = 6.6 x 10 10 cm). This late 
appearance of the high Lorentz factor component is consistent with the results by Zhang, 
Woosley, & MacFadyen (2003); Umeda et al. (2005). The maximum Lorentz factor in each 
model is in good agreement with the ones predicted by Eq. (2), (see, table 1 in which the 
maximum Lorentz factor (r maxsim ) achieved in the calculation and theoretically estimated 
value r max are presented.). Note that in model A300 the maximum Lorentz factor is expected 
to achieved after t — 10 s and is thus not shown in Fig. 7. When an outflow is injected to the 
computational domain, the radius of the outflow increases via adiabatic expansion. At this 
time the thermal energy is converted to the kinetic one and thus the Lorentz factor increases. 
This implies that the acceleration occurs due to the effective conversion of the thermal energy 
to kinetic energy during the propagation in the progenitor. As we show in Fig. 6, the Lorentz 
factor does not monotonically increase upward but decreases at the points where the density 
increases. These points correspond to the locations of oblique shocks. There exist not only 
oblique shocks but also rarefactions, in which the density decreases while the Lorentz factor 
increases. 

The pressure inside of the bow shock is very smooth except at the places where dis- 
continuities appear. These discontinuities follow each other, i.e. the separations between 
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them do not grow with time (see, Fig 8 which presents schematic figure of the positions 
of discontinuities). This feature is quite different from the shocks produced by a spherical 
explosion, such as supernova remnant, where the FS and RS separate with time. Most of 
the kinetic energy is converted to thermal energy at the RS. The hot gas between the RS 
and CD is exhausted to lateral direction and finally generates a back flow. 

Even after the breakout of the outflow from the progenitor surface we continued to 
calculate the propagation of the jet up to the distance of 6.5 x 10 10 cm. The outflow exhibits 
a significant expansion after the breakout. Fig. 9 shows rest mass density (upper) and 
Lorentz factor (lower) of model A50 at t — 4.75 s, which illustrates how the jet breaks 
out and expands to the ISM. Although the shocked progenitor gas and the cocoon which 
comprises the back flow are confined within a narrow zone by a strong bow shock before 
the breakout, this is no longer the case after the breakout. In fact, all the progenitor gas, 
including the back flow and the cocoon, begin to expand and forms a global outflow towards 
the ISM. Note, however, that still the high velocity component survives along the z-axis. As 
a result, high velocity component is surrounded by slow velocity and dense component. We 
will further follow jet propagation in near future. 



3.2. Expanding Outflows : Cases with Low r 

The outflows with smaller injection velocity behave very differently, compared with 
those with larger injection velocity presented above. In this section we discuss models G01, 
G10, and G50, in which the injection velocity is fixed to be v = 0.5c (r = 1.15). 

Figure 10 shows the contours of rest mass density and Lorentz factor of these three 
models at t — 7.5 s (G50), t = 8.0 s (G10), and t = 8.5 s (G01), respectively. In contrast 
with the previous cases the outflow shows an expanding structure like a fan. Let us remind 
that the injected outflow is initially parallel to the cylindrical axis. This expanding structure 
can be seen from the early phase of the evolution of the outflow. We also see that the RS 
separates from the FS with time in these cases. The high velocity region can be seen around 
the innermost region which ends at the reverse shock, forming a "disk" -like region with a 
large cross section (hereafter called as a disk), where the cross section stands for the area of 
the disk in the three dimensional space (see figure 5b and 8). Interestingly, the separation 
between the FS and RS increases with the time in the present cases, just as in supernovae 
remnant or pulsar wind, which makes a good contrast with the cases of the collimated jet. 
As the outflow expands with a large opening angle, the cross section of the CD increases. 
Since the cross section is proportional to square of the radius of the outflow, the flow with 
a large opening angle needs so large dragging power to proceed. This suppresses the back 
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flow formation at the head of the jet. As the cross section of the CD increases, the gas is 
gradually collected at the head of the flow because the FS sweeps up the progenitor gas. 
This works to enhance the expansion of the outflow to the lateral direction, thus an opening 
angle being increased. The half opening angle # S i m of the outflow at the break is 22° (G01), 
30° (G10), and 26° (G50), respectively. 

3.3. Intermediate Cases 

The continuous transition from the collimated jet to the expanding outflow takes place, 
as the Lorentz factor (and thus velocity) of the injected outflow decreases. This is nicely 
illustrated in Figure 11, which shows the density and T contours of models with a variety of 
To values. The specific internal energy of the injected outflow is fixed to be eo/c 2 = 0.1. For 
models A01, B01, C01, and D01 we find a collimated jet, the half opening angle of which 
is only a few degrees. The outflow starts to show an expanding tendency from model E01 
(v /c = 0.7, e /c 2 = 0.1), in which the half opening angle 9 sim ~ 12° at the time when the 
outflow breaks out, and finally model 101 exhibits a typical expanding-fan structure with 
9 ~ 32° at t — 10 s (see Table 1). The radius of the outflow gradually increases as the 
velocity of injected outflow decreases. The sideway expansion of the bow shock becomes 
significant when v is small, v /c < 0.7. 

The same transition also takes place in other series of models. We next fixed the internal 
energy to be eo/c 2 = 1 and eo/c 2 = 5. The opening angle in each model varies from a few 
degrees (model A10, and A50) to more than 30 degrees (models G50 and H10). The transition 
takes place in models D10 and D50, in which the specific internal energy eo/c 2 if fixed to 
1.0 and 5.0, respectively. We summarize these results in Fig. 3, in which circles, squares, 
and triangles indicate the cases with a half opening angle of the outflow at the time of the 
breakout or at the time t — 10 s being sim < 5°, 5 < 9 sim < 20°, and 9 sim > 20°, respectively. 
The RS separates from the CD and the FS in cases of the expanding outflow. The maximum 
Lorentz factor r max sim achieved in the calculation in every model is presented in Table 1. 
Those are good agreement with theoretical estimated r max . The numerical errors caused by 
the relativistic transverse velocity respect to the propagation direction described in the 2.2 
and A.l are only appreciable in the intermediate cases between collimated jet to expanding 
outflow. We estimate that the numerical errors produce an uncertainly of 10 percent in the 
Lorentz factor at which the transition happens. Our main conclusion is not affected by the 
numerical errors of this sort, however, since the numerical errors are negligible in outflow 
which keeps quite good collimation, as well as in non-relativistic outflows. 
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3.4. Dependence on the Energy Injection Point 

We do not know where and how the outflow forms and evolves in the progenitor at 
present. To see how sensitive the jet propagation properties are to the injection point, we 
also calculated model A50b, in which the lower boundary is set to be 2 x 10 7 cm from the 
center; that is, 1.8 x 10 8 cm closer than in other cases (2 x 10 8 cm). Other model parameters 
are the same as those of model A50. Since the progenitor has a very large density gradient 
in the radial direction at small radii, it is not trivial weather the outflow can drill this "wall" 
in this case. Note that Aloy et al. (2000) also set the inner boundary at 2 x 10 7 cm and 
deposited thermal energy around the boundary. 

Fig. 12 shows the results of model A50b at the time t — 6.5 s when the outflow breaks 
out. There is a difference in the early phase (< 3 s) of the dynamics compared with model 
A50. The higher density of progenitor gas prevents the outflow from proceeding to the 
direction of z axis. The outflow gradually drills the progenitor, and the bow shock spreads 
to lateral direction. Except for these subtle differences, the morphology and dynamics are 
very similar, once the outflow drills out the high density region. The outflow can keep the 
collimated structure during the propagation in the progenitor. A back flow appears in both 
cases. The outflow can finally break out the progenitor as shown in case A50. 



4. DISCUSSION 

4.1. Physics discriminating two types of outflow 

We observe a dramatic but smooth transition from collimated structure to expanding 
one in the outflow propagation by changing the Lorentz factor and specific internal energy 
of the injected outflow. What is the key physics which is responsible for the transition ? 
One of them is the pressure of the injected outflow. The pressure of the injected outflow 
in the models of the collimated jet is lower than that of expanding outflow. Additionally, 
we, here, point out that the appearance of the (internal) reconfinement shock could be a 
key, since the lateral expansion is suppressed by the presence of the reconfinement shock 
(Mizuta et al. 2004), thus leading to a formation of the collimated jet structure. In fact, 
high-pressure regions driven by a bow shock can keep the collimated structure. Models which 
have expanding structure have no or less internal structure in the reverse and side shocks. 
Then the injected flow is bi-forked. On the contrary, models which keep the collimated 
structure have a number of internal oblique shocks within the jet. The existence of a back 
flow also enhances the appearance of such shocks. Such shocks may allow the magnetic field 
generation in the jet (discussed later). 
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The internal shock model was introduced by Rees & Meszaros (1994 ; see also Kobayashi 
et al. 1997, and Spada et al. 2001) to explain very short time variation of GRBs. It 
is predicted that the outflow has up to a few hundred internal structures or "shells" and 
multiple two shell collisions occur at the distance of 10 13 ~ 15 cm from the fire ball which 
emits gamma-rays. 

A simple linear analysis by Urpin (2002) and Aloy et al. (2002) concludes that the 
timescale for a jet to propagate to the progenitor surface is long enough for perturbations to 
grow to form internal shocks and variability in GRBs. The instability is caused by the shear 
flow or perturbation in pressure, density, and/or velocity in lateral direction in the jet. The 
derived growth rate of the shear instability monotonically increases with the decrease of the 
wavelengths (Urpin 2002). Although with our current resolution it is impossible to resolve 
such fine structure, we can estimate the timescale of the instability, assuming the specific 
internal structure, say, that seen in Fig 6. Taking the wavelength to be a typical interval of 
the oblique shocks, ~ 10 9 cm, and bulk Lorentz factor to be ~ 5, we find that the timescale 
of the shortest mode in model A50 is on the order of 10~ 2 s. This growing timescale is 
sufficiently shorter than the dynamical one; that is, there is ample time for perturbations to 
grow to shocks, as was claimed by Aloy et al. (2002). 

In addition to the shear instability scenario, there exists another possible cause of the 
internal oblique shocks; that is the interaction between the jet and the back flow. As the 
path of the back flow is not straight, the boundary between main jet and the back flow is not 
smooth but dynamically perturbed. Since the dynamics is so complicated and is non-linear, 
it is hard to specify which mechanism, either the growth of the instability in the jet or the 
interaction between the jet and the back flow, is dominant for the appearance of the internal 
shocks. The detailed structure should depend on the resolution of the calculation. 

For this reason, we have calculated the same model as model A50 but with twice higher 
resolutions than that of model A50; that is model A50c. Fig. 13a shows the results of model 
A50c (see the bottom panel of Fig. 6 for comparison). The discontinuities still clearly appear 
in A50c, as well, but the number of discontinuity has increased than that of model A50. The 
cocoon has fine structures and vortices. The emergence of vortices enhances the mixing of 
the back flow and shocked progenitor gas. 

The appearance of the oblique shocks in the jet could be a key to the collimation. 
Whether such internal shocks appear or not strongly depends on how high pressure can be 
achieved by the presence of the bow shock. But it is difficult to predict whether the outflow 
expands or keeps collimated structure by a simple formula. What is clearly demonstrated 
through our simulations is that the large Lorentz factor and/or internal energy of the outflow 
can make a strong bow shock ahead of the outflow, leading to a good confinement of an 



-16- 



exploding hot outflow. 

Before closing this section, we wish to remark on the role of magnetic fields. The 
GRBs quite generally exhibit power-law and non-thermal spectra produced by synchrotron 
emission. Then, we need moderately strong magnetic fields within the jet, but it is not 
well understood how the magnetic field can be generated and grows in GRBs. A plausible 
mechanisms for creating magnetic fields is the Weibel instability which sets out when the 
velocity field is not isotropic and which can amplify magnetic fields quickly. Further, the 
particle acceleration to highly relativistic regime can take place in the shocks. Nishikawa et 
al. (2005) showed relativistic electrodynamics particle simulations of launching jet into the 
ambient medium. They observed the amplification of non-uniform and small-scale magnetic 
fields by the Weibel instability. Although they studied the shock at the head of the jet, 
magnetic fields can also be amplified in internal shocks within the jet. 

If the magnetic field grows up, it will inevitably affect the emissivity of pre-cursors 
and emissivity of gamma-rays of GRBs after breakout from the progenitor. If large-scale 
magnetic fields can be created, they should affect the dynamics of the outflow. We need 
further studies in this field. 



4.2. GRBs, X-ray Flashes, and SNe 

Finally we discuss the astrophysical implications of our results. We have shown differ- 
ent types of outflow dynamics which can arise even by the input of the same total energy 
power through the initial outflow. The input outflows may propagate, keeping a collimated 
structure in the progenitor, or showing an expanding structure. The change of the outflow 
shape with changes in e and r is gradual (see Fig. 3). Although we cannot directly observe 
the propagation of the outflow within the progenitor due to a large optical depth, different 
dynamics will produce observable effects after the outbreak of the outflow. Different outflow 
dynamics may account for a variety of phenomena. The outflow with significant collimation 
and high Lorentz factor (> 100) will produce GRBs. When the Lorentz factor in the colli- 
mated jet is a bit smaller, so is the peak energy of the emission. Such an outflow could be 
observed as XRFs. 

Some of our simulations show an expanding outflow even when T and/or eo are relatively 
small. Such outflows may be observed as less energetic explosions as XRFs, even if the 
viewing angle is relatively small. These cases may correspond to the high velocity but non- 
relativistic flow, such as SN2002ap, since the flow still has a directivity. 

There is another factor for explaining the differences between GRBs and XRFs; that is 
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the viewing angle to the jet. Recently the unified model has been proposed by Nakamura 
(2000) and by Yamazaki et al. (2002, 2004) who claim that the different phenomena are 
attributed to different viewing angles, although the bursts themselves are identical. If the 
viewing angle is less or near the jet opening angle, the object will be observed as a GRB. If 
not, it will be identified as an XRF. Zhang et al. (2004) derived the observational energy as 
a function of the viewing angle based on their simulations. 

The third but less unlikely possibility has been pointed out through the discovery of 
Soft Gamma-ray repeaters (SGR) 1806-20. Because of its proximity its spectral and light 
variations have been recorded in unprecedent details (Terasawa et al. 2005; Palmer et al. 

2005) . Such an event may be observed as a short gamma-ray bursts, if it explodes in nearby 
galaxies, although its spectrum is basically blackbody and does not agree with those of 
GRBs. 

Finally, we comment on the synthesis of heavy elements. The large cross section of the 
forward bow shock is advantageous for the synthesis of heavy elements, since then wider 
regions are available for the nucleosynthesis (MacFadyen and Woosley 1999; Nagataki et al. 
2003). We find that the outflows have an expanding structure in some cases. This expansion 
is reminiscent of the cases of aspherical explosions of supernova or hypernova as was shown in 
Nagataki et al. 2003 (1998, see also Nagataki 2000, Maeda & Nomoto 2003, Nagataki et al. 

2006) . A large amount of heavy elements is expected to be synthesized by such an aspherical 
expansion, since the effective cross section of the bow shock is very large. The difference of 
the expansion should affect the amount of the synthesized elements. The viewing angle is 
also an important parameter for this type of outflow and is expected to explain the observed 
features of SNe (Mazzali et al. 2005). 

Nagataki et al. (2003) calculated explosive nucleosynthesis by hydrodynamic calculations 
along the line of the collapsar model. It is also possible that the nucleosynthesis of heavy 
elements can occur in the accretion disk surrounding the black hole (Woosley, Eastman, 
& Schmidt 1999; MacFadyen and Woosley 1999; Fujimoto et al. 2001; Pruet et al. 2003). 
Nagataki et al. (2003) showed that 56 Ni is synthesized in the jet like outflow and concluded 
that the amount of 56 Ni becomes larger when the energy deposition occurs in a short time, 
since then the deposited energy is effectively converted to thermal energy by the strong 
shock. Such calculations can be a good diagnostic tool to compare with the observations, 
since the amount of synthesized elements can be estimated. 
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5. CONCLUSION 

We investigate the propagation and dynamics of the outflows in the progenitor in the 
context of collapsar model for a central engine of GRBs by means of hydrodynamical simu- 
lations. We assume a fixed power input of the initial outflow of E = 10 51 ergs s -1 and the 
radius of the injected outflow R = 7 x 10 7 cm, and follow the propagation of hot outflow 
for different values of e and r over the ranges of 1.05 < T < 5 (0.3 < v /c < 0.98) and 
0-1 < e /c 2 < 30. The net energy for 10-second injection satisfies the explosive energy of 
so-called hypernovae ~ 10 52 ergs. 

Our conclusions can be summarized as follows: 

1. The propagation dynamics of the outflow dramatically changes from the collimated 
structure to the expanding one as T decreases. If the Lorentz factor is high enough, 
say r > 3, the outflow can propagate throughout the progenitor, keeping a very 
collimated structure. The half opening angle is 6 S i m < 2° for To > 3. But the opening 
angle has weak dependence on eo, as well; we get 8 S i m < 3° even for smaller T but 
with small e /c 2 < 0.1. The maximum Lorentz factor is, on the other hand, sensitive 
to both of T and e ; roughly T max ~ r (l + e /c 2 ). 

2. In the relativistic, collimated flow, a back flow, which is anti-parallel to the main jet, 
appears. During the propagation in the progenitor, we can see some internal struc- 
tures caused by the instability grown by the shear flow in the jet or by the interaction 
between the jet and back flow in the collimated jets. Such oblique shocks can help 
the reconfinement of the jet. The maximum Lorentz factor of the jet follows a sim- 
ple formula derived from energy conservation relation. After the breakout the outflow 
expands into the interstellar space, although there still remains a high velocity compo- 
nent along the z-axis. Its half opening angle if a few degrees. This could be observed 
as GRBs. Another flow component which surrounds the central high velocity compo- 
nent can also be seen. It originates from the back flow during the propagation in the 
progenitor and shocked progenitor gas. 

3. When the Lorentz factor (or the initial velocity) of the outflow is not large, say, T < 1.4 
( v o/c ^ 0.7), the outflow no longer keeps the collimation and thus expands to the 
forward and lateral directions. Eventually the outflow breaks out like an aspherical 
supernova explosion. This is because with the small outflow velocity the bow shock 
is weak and cannot drive the progenitor gas to high enough pressure. As a result, 
the reconfinement shocks, which are necessary for the collimation, does not appear. 
Thus, the structure is relatively featureless in the outflow. As the cross section of the 
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reverse shock increases with time, the mass is collected at the head of the outflow. 
This enhances lateral expansion. 

4. High Lorentz factor (> 10)is needed to explain energetic phenomena, such as GRBs and 
XRFs, but the different initial internal energy affects the opening angle of the outflow 
for injected outflows of smaller Lorentz factor, i.e. slower velocity, thereby producing 
a marked difference in its observable. Rather low internal energy, eo/c 2 < 0.1, and 
relatively small Lorentz factor (< 5) leads to collimated non-relativistic jets, which will 
be observed as a failed GRB. High internal energy, eo/c 2 > 5, leads to un-collimated 
relativistic jets, which could be observed as XRFs. We can thus phenomenologically 
explain different types of explosions, GRBs XRFs, and failed GRB along the same 
line but with different values of e for slower injected velocity. However, a cause of 
producing a variety of e and r is still unknown. It should be investigated in future 
work. 
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A. ID and 2D Test Calculations 

Since the relativistic hydrodynamic code used in this study is updated version of the one 
used in Mizuta et al. (2004). In this appendix, we describe the numerical methods employed 
in the code, and show the results of ID and 2D numerical test problems that have been used 
by other groups to show the ability of our new code. 

Relativistic hydrodynamic equations can be written in conservative form as follows. 
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Cartesian coordinate: 



du | dfju) | dgju) _ Q 
<9t <9x 9y 



Cylindrical coordinate: 



du + 1 9(r/(u)) + 9gr(w) = 
dt r dr dz 

Spherical coordinate: 

du t 1 d(r 2 /(t*)) 1 <9(sinflc/(n)) _ 

where, u is conservative vector, / and g are flux vectors, s c and s s are source vectors, 
respectively. They are defined as follows, 

u = (pT, phT 2 v\ phT 2 v 2 , phT 2 -p- pT) T , (A4) 
f(u) = (pTv\ phY 2 v 1 v 1 + p, phT 2 v 1 v 2 ) phY 2 v 1 - pTv 1 ) 7 ', (A5) 
g(u) = (pTv 2 , phY 2 v 1 v 2 , phT 2 v 2 v 2 + p, phT 2 v 2 - pTv 2 ) T } (A6) 

s c = (0,^,0,0) T , (A7) 

2p + phT 2 v 2 v 2 phY 2 v V - cot 6p \ T 
0, , ,0 . (A8) 



In this appendix, we use the unit that speed of light is unity. We employ equation of state 
with constant adiabatic index {p— (7 — l)pe). The discretized formula for the version of the 
first order accuracy in time is as follows, 

u n+i = u n_ AtL(w n ). (A9) 

where subscripts indicate time steps, the time step is At, and a function of L{u) is defined 
as, for example in cylindrical coordinate case, 

L{u) = 

{ri+i/2j/(«)i+i/2j + n-yzjfiu)^^} + —{g{u) iij+1/2 - g(u) i:j _ 1/2 } + a c (w) i ^10) 



where subscripts stand for space. Time integration is carried out using TVD Runege-Kutta 
method (Shu & Osher 1989) to get higher order accuracy in time. For example, we use 
following formulae in the versions of the second and third order accuracy in time. At first 
common process is done, 

u {1) =u n + ML{u n ). (All) 
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is assigned to be u n+l for the version of first order accuracy in time as shown in Eq.A9. 
Following additional approaches are done for the version of the higher order accuracy in 
time. 

RK2 (2nd order accuracy) 

u n+1 = hu {1) + AtL(u {1) )) (A12) 

RK3 (3rd order accuracy) 

u (2) = ^ + + AtL( U M)) (A13) 

u n+l = \u n + \{u^ + AtL(u^)) (A14) 

The recovery from conservative values u to primitive variables, such as, p, p, v is done 
following Aloy et al. (1999). 

Our code is based on the Godunov type scheme which is applied in many relativis- 
ts hydrodynamic codes and has an advantage to catch the shocks clearly compared with 
other schemes. The piecewise constant, Monotone Upwind Scheme for Conservation Laws 
(MUSCL), and PPM are adopted for reconstruction to derive the cell surface states. The 
MUSCL reconstruction used in the code was described in Mizuta et al. (2004). We fol- 
lowed the PPM reconstruction as shown by Marti & Miiller (1996), see also original paper 
of the PPM by Collela & Woodward (1984) for non-relativistic hydrodynamic code. Those 
reconstruction can allow us to get 1st, 2nd, and 3rd order accuracy in space, respectively. 

Numerical fluxes, such as, f(u) and g(u) are derived from the Marquina's flux formula 
(Donat & Marquina 1996; Donat et al. 1998) instead of using an exact Riemann solver 
as Collela & Woodward (1984) and Marti & Miiller (1996) did. The expressions of the 
eigenvalues and left and right eigenvectors of Jacobian matrices, such as, df/du and dg/du 
which are necessary to calculate numerical fluxes are given in Donat et al. (1998). 

Recent development of not only relativistic hydrodynamic code but also relativistic 
magneto-hydrodynamic code allow us to simulate high energetic phenomena, such as rela- 
tivistic jets from AGN, quasars, and micro quasars, pulsar wind, and GRBs. The approach 
to relativistic hydrodynamic simulations by Godunov type scheme was employed in 90 's. 
Marti & Miiller (1996) developed the PPM code which uses an exact Riemann solver of 
relativistic hydrodynamics. Eulderink & Mellema (1995) derived the Roe average for the 
general relativistic hydrodynamic equations and extended the original Roe scheme, which is 
for non-relativistic hydrodynamics, to general relativistic hydrodynamics. Recently a num- 
ber of codes have been developed based on Godunov-type scheme which uses an approximate 
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Riemann solver to derive the numerical fluxes for not only relativistic hydrodynamics (Dun- 
can & Hughes (1994); Font et al. (1994); Donat et al. (1998); Del Zanna & Bucciantini 
(2002); Lucas-Serrano et al. (2004); Mignone k Bodo (2005); Zhang & MacFadyen (2005); 
Rahman & Moore (2005)), but also relativistic magneto- hydrodynamics (Komissarov (1999); 
Balsara (2001); Gammie et al. (2003); Leismann et al. (2005); Shibata & Sekiguchi (2005); 
Anton et al. (2006)), see also useful review by Marti Miiller (2003). The codes which do 
not use an exact or approximate Riemann solver have also been developed by van Putten 
(1991); Koide (2003); De Villiers & Hawley (2003). 

The results of the following all test calculations are done by the version of the PPM 
reconstruction and a third order accuracy version by TVD Runge-Kutta for time integration. 
The same version is used for our main results, i.e, 2D outflow propagation in the collapsar. 
The detailed analysis of test results, such as the error by different accuracy and resolution 
will be presented in the paper in prep. 

A.l. ID shock tube problem with no transverse velocity 

The shock tube problem is one of the standard test problems for not only the compressive 
non-relativistic hydrodynamic codes but also relativistic hydrodynamic codes, since we can 
test the features of shocks, rarefactions, and contact discontinues by this problem and the 
analytic solutions are available via iterative calculations. The problem is a kind of Riemann 
problem. Two half finite constant states are assumed as an initial condition (t = 0). We have 
tested two standard sets of initial parameters which were tested by the most of relativistic 
hydrodynamic codes. These cases do not include any transverse velocity. Those are shock 
tube problems A and B, and detailed parameters, geometry and total grid points are : 

• Shock tube A (plane) : 400 uniform grid points 

Left state (0 < x < 0.5); pi = 10, pl = 13.3, v xL = 0,v VL = 0,7l = 5/3 
Right state (0.5 < x < 1); p R = l,p R = 1 x lO" 6 ,^ = 0,v VL = 0,>y R = 5/3 

• Shock tube B (plane) : 400 uniform grid points 

Left state (0 < x < 0.5); pl — 1,Pl — 1000, v xL = 0,v VL = 0,7^ = 5/3 
Right state (0.5 < x < 1); p R = l,p R = 0.01, v xR = 0,v yR = 0,^ R = 5/3 

Figure 14 shows numerical and analytic solutions of these problems at t — 0.5 (shock 
tube A) and t = 0.35 (shock tube B). The profiles of the density, pressure, and 3- velocity 
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are presented. The numerical results of both problems are quite good agreement with the 
analytic solution, especially rarefactions. In case A, shock front is captured in good accuracy 
by several grid points. The contact discontinuity is a little diffusive. In case B, because of 
large pressure jump of initial condition, the jump of the density is about 10 times at the 
shock and very narrow shocked region appears. It is hard to resolve such narrow region by 
tested resolution. But our results are within the similar level of the results presented by 
other groups. 



A. 2. ID shock tube problem with transverse velocity 

If there is non-zero transverse velocity for ID shock tube problem, the results changes 
because of the dependence of the Lorentz factor on the absolute value of the velocity, as 
discussed by Pons et al. (2000); Rezzolla & Zanotti (2001). Recently it has been reported 
that it is numerically hard to resolve such problems with as same level of grid points as the 
cases without transverse velocities Mignone & Bodo (2005); Zhang & MacFadyen (2005). 
We followed the conditions by Mignone & Bodo (2005) in which the effect of some different 
initial transverse velocities cases was presented. The initial conditions are: 

• Shock tube C (plane) uniform 400 grid points 

Left state (0 < x < 0.5); p L = 1.0, p L = 1000.0, v xL = 0,7 L = 5/3 
Right state (0.5 < x < 1); p R = 1.0, p R = 1 x lO -2 ,-^ = 0,7^ = 5/3 

We have done parametric study by changing v VL from to 0.99 and v yR from to 0.99 as 
Mignone & Bodo (2005) did. Figure 15 shows the results of these tests. The profiles of the 
density, pressure and x-component of the 3-velocity at t = 4.0 are presented. The top-left 
panel (case (v VL , v yR ) = (0,0)) corresponds to the case Shock tube B presented above, but the 
result at different time is shown. The rarefactions are resolved in any cases as shown the 
cases without transverse velocity (shock tube A and B). As v VL increases, both positional and 
absolute value errors at around the discontinuities, such as, shocks and contact discontinuities 
increase. 

We also done resolution study on the set of (v VL , t> yi? ) = (0.9,0.9) from uniform 400 zones 
to 12800 zones. This has been done by Zhang & MacFadyen (2005) but they used their 
Adaptive Mesh Refinement (AMR) version for this. Figure 16 shows numerical and analytic 
solutions of this problem. The rarefaction is resolved in good accuracy by not only higher 
resolution calculation but also lower one. On the contrary, both the shock and the contact 
discontinuities are not resolved in both position and absolute value by lower resolution cal- 
culations. Higher resolution calculations can allow us to resolve shock front within a few 
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percents error in position. On the contrary, there still remains about several percents error 
in v y at the contact discontinuity. One of solutions to avoid such errors and to get good 
accuracy is to employ the AMR method which uses locally higher resolution at the area 
where discontinuities, such as, shocks and contact discontinuities exist to save the CPU time 
and memory requirement. The AMR was employed by Zhang & MacFadyen (2005) is one 
of the ways to reduce the error seen in this problem. 



A. 3. ID Reflection Shock Problem 

Reflection shock problem is suitable to test strong shocks. At t — 0, an uniform, high 
Mach number, and cold flow reflects at the boundary (x = or r = 0) which is a wall for the 
Cartesian coordinate case, a cylindrical axis for the cylindrical coordinate case, and center 
of the spherical coordinate for the spherical coordinate case. A strong shock appears and 
proceeds to the cold gas. The fluid is heated by this shock and becomes at rest. The shock 
front satisfies Rankine-Hugoniot relation for relativistic hydrodynamics (Taub 1948). An 
analytic solution is available (Johnson & McKee 1971), if the gas of the initial flow is cold 
(Po ~ 0). 

The analytic solution at t — T is, 

7 + 1 7 



7 — 1 7 — 1 



(r - 1) Po, e = To - 1, v = 0, for x < V S T or r < V S T (A15) 



x 



a 



p = po ( 1 + — ) , e ~ 0, v = v , for x > V S T or r > V S T, (A16) 



where V s is shock velocity, 



(V - 1 \ 1/2 

v -=(j&) < A17 > 

The rest mass density, specific internal energy, velocity, and Lorentz factor of the inflow are 
Po, eo 5 v o, and To, respectively. The geometry is plane (a = 0), cylindrical (a = 1), and 
spherical (a = 2), respectively. The expression of rest mass density jump at the shock front 
is divided in two parts, namely a maximum density compression ratio in non relativistic 
hydrodynamics ((7 + l)/(7 — 1)) and including a Lorentz factor. 

• REP 

p = 1.0, e = l0~ A ,v = -0.999(r = 22), 7 = 4/3, Plane 

• REC 

p = 1.0, e = 10" 4 ,t; = -0.999(r = 22), 7 = 4/3, Cylindrical 



-25 - 



• RES 

p = i.o, eo = W~ 4 ,v = -0.999(r = 22), 7 = 4/3, Spherical 

Figure 17 (a)-(c) show numerical and analytic solutions of this problem at t — 1.57. The 
calculations are done uniform 400 zones in the Cartesian (a), cylindrical (b), and spherical 
coordinate (c), respectively. Since we have taken care of the treatment of the boundary for 
the cylindrical and spherical coordinate cases, the error around the boundary (singular) is 
reduced compared with the results shown in Mizuta et al. (2004). The error at and round 
the wall boundary is studied by Noh (1987) in detail. He discussed not only case in plane 
geometry but also cylindrical and spherical geometry. Our numerical error at the wall is 0.8 
% (Cartesian), 1.4% (cyrindrical), and 8.3% (spehrical). 

Figure 17 (d) shows numerical and analytic solutions of reflection shock problem by 
Cartesian coordinate, but different initial velocity, from —0.9 up to —0.999999999, corre- 
sponding Lorentz factor is from 2.27 to 22361. Uniform 400 zones are used. The other 
conditions, such as, density and pressure is as same as the problem (REP). In all cases, the 
strong shock front is captured with a few grid points. All results are at t — 2. The jump 
condition at the shock and location of the shock represents analytic solution. 



Two dimensional shock tube problem is done to confirm the shock dynamics in the 
multidimensional case. This problem includes the interactions of shocks, rarefactions, contact 
discontinuities. Initially a square computational domain is prepared in x-y plane and divided 
into four quarter boxes (see left panel in Fig. 18). Initial conditions in each box are : 



Adiabatic index is assumed to be 7 = 5/3. This set of initial condition is as same as used by 
Del Zanna & Bucciantini (2002); Zhang & MacFadyen (2005). This is also similar condition 
done by Lucas-Serrano et al. (2004). At first two shocks appear from the high pressure 
regions A and D and they proceed into the region B. These shocks collide each other in the 
region B. As a result, a high pressure gas break into the region C. We use 400 x 400 uniform 
grid points in a square computational box. The boundary conditions at all boundary are 
open ones. 



A. 4. 2D shock tube problem 



(p,V X ,Vy,p) 
(p,V x ,Vy,p) 
(p,V x ,Vy,p) 
(p,V x ,Vy,p) 



(0.10,0.00,0.00,0.01) : 0.5 < x < 1, 0.5 < y < 1 (region A), 
(0.10, 0.99, 0.00, 1.00) : < x < 0.5, 0.5 < y < 1 (region B), 
(0.50, 0.00, 0.00, 1.00) : < x < 0.5, < y < 0.5 (region C), 
(0.10, 0.00, 0.99, 1.00) : 0.5 < x < 1, < y < 0.5 (region D). 
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Figure 18 shows 30 levels of iso-surface of the logarithm of rest mass density at t — 0.4. 
Two curve of shocks in the region B are well resolved as shown by other groups. It should 
be noted that sound waves appear at the contact discontinuities between the region A and 
C, and the region C and D can be seen in the density contour plot in the region C. Such 
wave can also be seen in the results Del Zanna & Bucciantini (2002); Lucas-Serrano et al. 
(2004); Zhang & MacFadyen (2005). 

A. 5. 2D Double Mach Reflection Problem 

2D double Mach reflection problem was introduced by Woodward & Colella (1984) for 
non-relativistic hydrodynamics and recently has been applied to relativistic hydrodynamics 
by Zhang & MacFadyen (2005). An initial uniform shock collides with a reflective wall, 
forming another shock. These two shocks interact each other. A kind of jet appears along 
the wall. The solution is self similar. 

We have tested this problem, using as same parameters as Zhang & MacFadyen (2005) 
did. The density and pressure of the unchecked gas is 1.4 and 0.0025, respectively. The 
classical shock Mach number, M c = vs/cs, where vs is shock velocity and cs is sound 
velocity of the unchecked gas, respectively, is 10. The adiabatic index is 7 = 1.4. The shock 
velocity in this case is ~ 0.4984 accordingly. The state of initial shocked gas can be derived by 
Rankine-Hugoniot relation of relativistic hydrodynamics. The density, pressure, and velocity 
of shocked gas are ~ 8.564, ~ 0.3808, and ~ 0.4247, respectively. The computational domain 
is x — y plane and 4x1 rectangular which includes uniform 512 x 128 zones. The shock 
front makes an angle of 60 degrees with x axis and the rim of the shock is at x — 1/6, y = 
initially. The boundary conditions at x = 0, 0<£<l/6aty = 0, and < x < V s sin 60° t 
at y — 1 are inflow of the initial shocked gas. The reflective boundary condition is employed 
at 1/6 < x < 4,y = 0. The boundary condition at other boundaries are open condition. 

Figure 19 shows 30 levels of iso-contours of the rest mass density at t — 4.0. The global 
feature of the shocks is represented as shown by Zhang & MacFadyen (2005). Although the 
resolution of presented calculation is not so high, the Kelvin-Helmholtz instability which 
violates the self similarity of the solution can be seen at around x = 2.5, y — 0. 

A. 6. 2D Emery Step Problem 

The problem of a wind tunnel containing a step was proposed by Emery (1968) has been 
tested (2D Emery problem). He tried to compare the results by different schemes. A step 
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boundary is introduced in the rectangular computational box (x x y = 3 x 1). This step is 
located at 0.6 from left boundary of the computational box and the hight is 0.2 and continues 
to the right boundary. The left side boundary is an inflow as the same quantities. Initially 
computational box is filled with a uniform supersonic flow (v x = 0.999, v y — 0, p — 1.4, and 
Newtonian Mach number is 3). The outflow (zero gradient) boundary condition is employed 
at right side of the boundary (x — 3). The reflective boundary conditions are employed at 
the step boundaries, y — 1, and < x < 0.6 at y = 0. 120 x 40 uniform grid points are 
spaced. The adiabatic index for equation of state is 7 = 1.4. The corner of the step becomes 
singular. The corrections around the corner of the step have been included, see Appendix 
in Donat et al. (1998). 

Figure 20 shows 30 levels of iso-surface of the logarithm of rest mass density. The global 
features, such as shocks and rarefactions, are reproduced as shown in Zhang & MacFadyen 
(2005). 

A. 7. Spherical Blast Wave 

At last, we show the results of spherical blast wave, using cylindrical coordinate. This 
problem is done to see how spherical symmetry is kept using cylindrical coordinate. The 
computational box is0<z<l,0<r<l with equally spaced grid points (320 x 320). 
Initially, a high pressure gas (p = 1000) is put in the region, \J z 2 + r 2 < 0.4. The pressure 
in the outside is p — 1. At t = all gas is at rest and the density is uniform p — 1. The 
adiabatic index for equation of state is 7 = 5/3. The high pressure gas expands and a 
shock proceeds into cold gas. The boundary conditions are reflective one at both axises and 
zero gradient outflow one at the other boundaries. The same problem has been done using 
spherical coordinate with 3200 uniform zones for comparison. 

Figure 21 shows one dimensional rest mass density profile of along both z (top) and r 
(bottom) axises at t — 0.4. The solid lines are the results of the same problem by spherical 
coordinate with 3200 uniform zones in radius. The both results along z and r axises are in 
good agreements with spherical one, although the velocity profile along z axis has a bump 
around the head of the wave x ~ 0.75. 
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Table 1: Calculated models. The basic model parameters are the Lorentz factor (r ) and the 
specific internal energy (e ). Also shown are the ratio (/ ) of the kinetic energy to the total 
energy (where rest mass energy is excluded) and the estimated maximum Lorentz factor 
(r m ax)- The achieved maximum Lorentz factor (r maXiSim ) and half opening angle (# s im) from 
the calculations are also listed. 



model 


eo/c 


J- o [vo/c) 


Jo 


f - 1 
Jo 


1 max 


J- max,sim 


a 

"sim 


A01 


0.1 




0.86 


1.2 


5.5 


5.6 


1.2 


A10 


1.0 


5 (0.98) 


0.38 


2.7 


10 


11 


1.2 


A50 


5.0 




0.11 


9.3 


30 


31 


1.7 


A300 


30 




0.020 


51 


155 


104 


1.0 


B01 


0.1 




0.85 


1.2 


A A 

4.4 


A V 

4.5 


1.6 


BIO 


1.0 


4 (0.97) 


0.36 


2.8 


8.0 


9.0 


1.7 


B50 


5.0 




0.10 


9.8 


24 


26 


1.3 


C01 


0.1 




0.84 


1.2 


3.3 


O A 

3.4 


1.5 


CIO 


1.0 


3 (0.94) 


0.34 


2.9 


6.0 


6.7 


1.8 


C50 


5.0 




0.093 


11 


18 


on 

20 


1.9 


1JU 1 


n 1 

U. 1 




U.oU 


1 ^ 

l.O 


9 9 
_._ 


9 ^ 


9 Q 


D10 


1.0 


2 (0.87) 


0.29 


3.5 


4.0 


4.5 


15 


D50 


5.0 




0.074 


14 


12 


14 


12 


E01 


0.1 




0.71 


1.4 


1.5 


1.6 


12 


E10 


1.0 


1.4 (0.7) 


0.20 


5.1 


2.8 


3.1 


20 


E50 


5.0 




0.0470 


21 


8.4 


9.4 


19 


F01 


0.1 




0.64 


1.6 


1.4 


1.4 


19 


F10 


1.0 


1.25 (0.6) 


0.15 


6.6 


2.5 


2.8 


28 


F50 


5.0 




0.034 


29 


7.5 


8.3 


26 


G01 


0.1 




0.55 


1.8 


1.3 


1.4 


22 


G10 


1.0 


1.15 (0.5) 


0.11 


9.1 


2.3 


2.6 


30 


G50 


5.0 




0.024 


41 


6.9 


7.7 


26 


HOI 


0.1 


1.1 (0.4) 


0.44 


2.3 


1.2 


1.3 


27 


H10 


1.0 




0.073 


14 


2.2 


2.5 


29 


101 


0.1 


1.05 (0.3) 


0.31 


3.2 


1.2 


1.2 


32 



-33- 




1e+10 2e+10 3e+10 4e+10 5e+10 
radius cm 



Fig. 1. — The adopted radial mass density profile of the progenitor from the center (left) of 
the core to the surface (right). We adopt the profile at r > 2 x 10 8 cm as our initial mass 
profile, assuming spherical symmetry. This is equivalent to assuming that a mass of about 
two solar masses has collapsed. 
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Fig. 2. — A schematic figure of the progenitor and computational box throughout the present 
study. The radius of the injected outflow (i? ) is 7x 10 7 cm. The distance between the center 
of the progenitor and the lower-left corner of the computational box (z\ ow ) is 2 x 10 8 cm or 
2 x 10 7 cm. 
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Fig. 3. — The contours of the log-scaled mass density (in g cm -3 ) on the (To, eo/c 2 )-plane. 
The total energy flux is fixed to be Eq = 10 51 ergs s -1 and the radius of the injected outflow 
is i?o — 7 x 10 7 cm. The circles, squares, and triangles correspond to models, in which 
the half opening angle of the outflow at the time of breakout or at t = 10 s is 6 s[m < 5°, 
5° < sim < 20°, and # sim > 20°, respectively. 



Fig. 4. — The contours of the rest mass density (upper) and the Lorentz factor (lower) of 
models, from top to bottom, A300, A50,, A10, and A01, in which the Lorentz factor of 
injected outflow r is fixed to be 5, at the time of the breakout; t = 3.50 s (A300), t = 3.50 
s (A50), t = 3.00 s (A10), and t = 2.75 s (A01), respectively. All models keep collimated 
structure. A back flow, which runs from the head of the jet in the anti-parallel direction, 
appears in each model. 
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(a) Collimated Jet 

Bow Shock 




(b) Expanding Outflow 




(disk in 3D) 

Fig. 5. — Schematic figures of the collimated flows (a) and the expanding outflow (b). Only 
the main flow (outflow and back flow) is displayed here. The back flow is usually meandering. 
The half opening angle 8 S i m which is measured when the jet breaks out of the progenitor is 
presented. 
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z [cm] 

Fig. 6. — One dimensional profiles of density (solid line) and Lorentz factor (dashed line) 
along the z axis in the early phase of each simulation. Several discontinuities can be seen. 
Those correspond to the oblique shocks in the jet. For the hotter outflow, Lorentz factor 
increases at the oblique shocks. 
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Fig. 7. — One dimensional profiles of density (top) and Lorentz factor (bottom) along the z 
axis of model A300 at different time (t = 2, 4, 6, 8 and 10). Maximum Lorentz factor is a 
few tens during propagation in the progenitor and increases larger than 100 after breakout 
from the progenitor. 




Fig. 8. — A figure of the positions of froward shock, contact discontinuity, and reverses shock 
for cases of collimated jet (solid line) and expanding outflow (dashed line). 
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Fig. 9. — Same as Fig. 4 but for model A50 at t — 4.75 s. All the flow components, including 
the back flow during the propagation in the progenitor, become outflow. The region with 
high velocity gas remains along the z-axis. 
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Fig. 10. — Same as Fig. 4 but for models G50 (top), G10 (middle), and G01 (bottom), in 
which the velocity of injected outflow v is fixed to be 0.5c, at t — 7.5 s (G50), t = 8.0 s 
(G10), and t = 8.5 s (G01), respectively. The outflow has an expanding structure like a fan. 
As time goes on, the distance between the reverse shock (RS) and the forward shock (FS) 
increases, just as in a supernova remnant. 



Fig. 11— Same as Fig. 4 but for models A01, B01, C01, D01, E01, F01, G01 HOI, and 101. The transition in the 
outflow dynamics from the collimated structure to the expanding structure is evident. 
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Fig. 12. — Same as rest Fig. 4 but for model A50b in which the inner boundary is set to be 
at 2 x 10 7 cm from the center. In the beginning of the evolution it takes slightly a longer 
time to drill high density region. After that the dynamics and morphology are very similar 
to those of the regular case (A50). 
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Fig. 13. — (a) One dimensional rest mass density (solid line) and Lorentz factor (dashed 
line) profiles along the z-axis of model A50c, which is twice higher resolution version of 
model A50. Fine structure can be seen (see the result shown in the bottom of the Fig. 6 for 
comparison.), (b) Same as Fig. 4 but for model A50c. Complex structures can been seen in 
the cocoon caused by the mixing of back flow and shocked progenitor gas. 




Fig. 14. — Numerical (points) and analytic (solid lines) solutions of ID shock tube problem 
without transverse velocity cases: shock tube A (left) at t = 0.5 and shock tube B (right) at 
t = 0.35. Density, pressure and velocity profiles are shown. Uniform 400 zones are used for 
the calculation. 
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Fig. 15. — Numerical (points) and analytic (solid lines) solutions of shock tube problem with transverse velocity case 
with uniform 400 zones (shock tube C) at t — 0.4 are presented. From left to right, v yR = 0, 0.9, 0.99 and from top to 
bottom v yL = 0, 0.9, 0.99. Density, pressure, x-component of velocity are shown. 
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Fig. 16. — Numerical (points) and analytic (solids) solutions of shock tube problem with 
transverse velocity case (shock tube C, (v yL ,v yR )— (0.9,0.9)): Different resolution of the 
calculations are presented, from top to bottom, the number of grid points are 400, 800, 1600, 
3200, 6400, 12800, respectively. Density, pressure and x- and y-velocity components are 
shown. 
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Fig. 17. — Numerical (points) and analytic (solid lines) solutions of reflection shock problem 
with 400 uniform zones. The calculations are done in (a) plane geometry, (b) cylindrical 
geometry, (c) spherical geometry (t> = —0.99). The results at t — 1.57 are presented. The 
results of rest mass density for different initial velocity from v = —0.9 to v = —0.999999999 
at t — 2.0 are shown in (d). 
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Fig. 18. — Schematic figure of the initial condition of 2D shock tube problem (left) and 
numerical results (right), uniform 400 x 400 zones are used. 30 levels of iso-surface of the 
logarithm of rest mass density at t — 0.4 is shown. 
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Fig. 19. — 30 levels of iso-surface of the rest mass density of double Mach shock at t — 4.0. 
512 x 128 uniform zones and adiabatic index 7 = 1.4 for equation of state are used. 




Fig. 20. — Emery step at t — 4.0 with 120 x 40 uniform zones. 30 levels of iso-surface of the 
logarithm of rest mass density is shown. 
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Fig. 21. — Spherical blast wave problem by cylindrical coordinate with uniform 320 x 320 
(r x z) zones. Rest mass density, pressure, velocity along the axis at t — 0.4 are presented. 
Top: along z axis. Bottom: along r axis. Solid lines are results of the same problem but 
done by spherical coordinate with uniform 3200 zones. 



